## NULL
Show revision, parameter and session details
## Git revision: 4bc7e84
## git diff:
##  M processed_data/DDE_samples_as_replicates.tsv.gz
##  M processed_data/DE_per_cell_type_conditioned_on_bulk.tsv.gz
##  M processed_data/DE_samples_as_replicates.tsv.gz
##  M processed_data/aggregated_expression.tsv.gz
##  M processed_data/dde_one_vs_sum_of_rest_dt.tsv.gz
##  M processed_data/meta_and_reductions.tsv.gz
##  M scRNAseq_figures_using_only_processed_data.Rmd
##  M scRNAseq_figures_using_raw_data.Rmd
## ?? analysis_pipeline/
## ?? results/
## ?? scRNAseq_figures_using_only_processed_data.html
## ?? scRNAseq_figures_using_only_processed_data.md
## ?? scRNAseq_figures_using_only_processed_data_files/
## ?? scRNAseq_figures_using_raw_data.html
## ?? scRNAseq_figures_using_raw_data.md
## ?? scRNAseq_figures_using_raw_data_files/
## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 4.1.0 (2021-05-18)
##  os       CentOS Linux 7 (Core)       
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       Europe/Berlin               
##  date     2021-12-01                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
##  package          * version    date       lib source             
##  abind              1.4-5      2016-07-21 [1] standard (@1.4-5)  
##  AnnotationDbi      1.56.1     2021-10-29 [1] standard (@1.56.1) 
##  assertthat         0.2.1      2019-03-21 [1] CRAN (R 4.1.0)     
##  backports          1.3.0      2021-10-27 [1] standard (@1.3.0)  
##  Biobase            2.54.0     2021-10-26 [1] standard (@2.54.0) 
##  BiocFileCache      2.2.0      2021-10-26 [1] standard (@2.2.0)  
##  BiocGenerics       0.40.0     2021-10-26 [1] standard (@0.40.0) 
##  biomaRt            2.50.0     2021-10-26 [1] standard (@2.50.0) 
##  Biostrings         2.62.0     2021-10-26 [1] standard (@2.62.0) 
##  bit                4.0.4      2020-08-04 [1] standard (@4.0.4)  
##  bit64              4.0.5      2020-08-30 [1] standard (@4.0.5)  
##  bitops             1.0-7      2021-04-24 [1] standard (@1.0-7)  
##  blob               1.2.2      2021-07-23 [1] standard (@1.2.2)  
##  broom              0.7.10     2021-10-31 [1] standard (@0.7.10) 
##  bslib              0.3.1      2021-10-06 [1] standard (@0.3.1)  
##  cachem             1.0.6      2021-08-19 [1] standard (@1.0.6)  
##  cellranger         1.1.0      2016-07-27 [1] standard (@1.1.0)  
##  cli                3.1.0      2021-10-27 [1] standard (@3.1.0)  
##  cluster            2.1.2      2021-04-17 [2] CRAN (R 4.1.0)     
##  codetools          0.2-18     2020-11-04 [2] CRAN (R 4.1.0)     
##  colorspace       * 2.0-2      2021-06-24 [1] standard (@2.0-2)  
##  cowplot            1.1.1      2020-12-30 [1] standard (@1.1.1)  
##  crayon             1.4.2      2021-10-29 [1] standard (@1.4.2)  
##  curl               4.3.2      2021-06-23 [1] standard (@4.3.2)  
##  data.table       * 1.14.2     2021-09-27 [1] standard (@1.14.2) 
##  DBI                1.1.1      2021-01-15 [1] standard (@1.1.1)  
##  dbplyr             2.1.1      2021-04-06 [1] standard (@2.1.1)  
##  deldir             1.0-6      2021-10-23 [1] standard (@1.0-6)  
##  digest             0.6.28     2021-09-23 [1] standard (@0.6.28) 
##  dplyr            * 1.0.7      2021-06-18 [1] standard (@1.0.7)  
##  ellipsis           0.3.2      2021-04-29 [1] standard (@0.3.2)  
##  evaluate           0.14       2019-05-28 [1] standard (@0.14)   
##  fansi              0.5.0      2021-05-25 [1] standard (@0.5.0)  
##  farver             2.1.0      2021-02-28 [1] standard (@2.1.0)  
##  fastmap            1.1.0      2021-01-25 [1] standard (@1.1.0)  
##  filelock           1.0.2      2018-10-05 [1] CRAN (R 4.1.0)     
##  fitdistrplus       1.1-6      2021-09-28 [1] standard (@1.1-6)  
##  forcats          * 0.5.1      2021-01-27 [1] standard (@0.5.1)  
##  fs                 1.5.0      2020-07-31 [1] standard (@1.5.0)  
##  future             1.23.0     2021-10-31 [1] standard (@1.23.0) 
##  future.apply       1.8.1      2021-08-10 [1] standard (@1.8.1)  
##  gamlss             5.3-4      2021-03-31 [1] standard (@5.3-4)  
##  gamlss.data        6.0-1      2021-03-18 [1] standard (@6.0-1)  
##  gamlss.dist        5.3-2      2021-03-09 [1] standard (@5.3-2)  
##  gamlss.tr          5.1-7      2020-07-13 [1] standard (@5.1-7)  
##  generics           0.1.1      2021-10-25 [1] standard (@0.1.1)  
##  GenomeInfoDb       1.30.0     2021-10-26 [1] standard (@1.30.0) 
##  GenomeInfoDbData   1.2.7      2021-11-02 [1] standard (@1.2.7)  
##  ggforce            0.3.3      2021-03-05 [1] standard (@0.3.3)  
##  ggplot2          * 3.3.5      2021-06-25 [1] standard (@3.3.5)  
##  ggrepel            0.9.1      2021-01-15 [1] standard (@0.9.1)  
##  ggridges           0.5.3      2021-01-08 [1] standard (@0.5.3)  
##  globals            0.14.0     2020-11-22 [1] standard (@0.14.0) 
##  glue               1.5.0      2021-11-07 [1] standard (@1.5.0)  
##  goftest            1.2-3      2021-10-07 [1] standard (@1.2-3)  
##  gridExtra          2.3        2017-09-09 [1] standard (@2.3)    
##  gtable             0.3.0      2019-03-25 [1] standard (@0.3.0)  
##  haven              2.4.3      2021-08-04 [1] standard (@2.4.3)  
##  hms                1.1.1      2021-09-26 [1] standard (@1.1.1)  
##  htmltools          0.5.2      2021-08-25 [1] standard (@0.5.2)  
##  htmlwidgets        1.5.4      2021-09-08 [1] standard (@1.5.4)  
##  httpuv             1.6.3      2021-09-09 [1] standard (@1.6.3)  
##  httr               1.4.2      2020-07-20 [1] standard (@1.4.2)  
##  ica                1.0-2      2018-05-24 [1] standard (@1.0-2)  
##  igraph             1.2.8      2021-11-07 [1] standard (@1.2.8)  
##  IRanges            2.28.0     2021-10-26 [1] standard (@2.28.0) 
##  irlba              2.3.3      2019-02-05 [1] standard (@2.3.3)  
##  jquerylib          0.1.4      2021-04-26 [1] standard (@0.1.4)  
##  jsonlite           1.7.2      2020-12-09 [1] standard (@1.7.2)  
##  KEGGREST           1.34.0     2021-10-26 [1] standard (@1.34.0) 
##  KernSmooth         2.23-20    2021-05-03 [2] CRAN (R 4.1.0)     
##  knitr              1.36       2021-09-29 [1] standard (@1.36)   
##  later              1.3.0      2021-08-18 [1] standard (@1.3.0)  
##  lattice            0.20-45    2021-09-22 [1] standard (@0.20-45)
##  lazyeval           0.2.2      2019-03-15 [1] standard (@0.2.2)  
##  leiden             0.3.9      2021-07-27 [1] CRAN (R 4.1.0)     
##  lifecycle          1.0.1      2021-09-24 [1] standard (@1.0.1)  
##  listenv            0.8.0      2019-12-05 [1] standard (@0.8.0)  
##  lmtest             0.9-39     2021-11-07 [1] standard (@0.9-39) 
##  lubridate          1.8.0      2021-10-07 [1] standard (@1.8.0)  
##  magrittr           2.0.1      2020-11-17 [1] standard (@2.0.1)  
##  MASS               7.3-54     2021-05-03 [2] CRAN (R 4.1.0)     
##  Matrix             1.3-4      2021-06-01 [2] CRAN (R 4.1.0)     
##  matrixStats        0.61.0     2021-09-17 [1] standard (@0.61.0) 
##  memoise            2.0.0      2021-01-26 [1] standard (@2.0.0)  
##  mgcv               1.8-38     2021-10-06 [1] standard (@1.8-38) 
##  mime               0.12       2021-09-28 [1] standard (@0.12)   
##  miniUI             0.1.1.1    2018-05-18 [1] standard (@0.1.1.1)
##  modelr             0.1.8      2020-05-19 [1] standard (@0.1.8)  
##  munsell            0.5.0      2018-06-12 [1] standard (@0.5.0)  
##  nlme               3.1-153    2021-09-07 [1] standard (@3.1-153)
##  parallelly         1.28.1     2021-09-09 [1] standard (@1.28.1) 
##  patchwork        * 1.1.1      2020-12-17 [1] standard (@1.1.1)  
##  pbapply            1.5-0      2021-09-16 [1] standard (@1.5-0)  
##  pillar             1.6.4      2021-10-18 [1] standard (@1.6.4)  
##  pkgconfig          2.0.3      2019-09-22 [1] standard (@2.0.3)  
##  plotly             4.10.0     2021-10-09 [1] standard (@4.10.0) 
##  plyr               1.8.6      2020-03-03 [1] standard (@1.8.6)  
##  png                0.1-7      2013-12-03 [1] standard (@0.1-7)  
##  polyclip           1.10-0     2019-03-14 [1] standard (@1.10-0) 
##  prettyunits        1.1.1      2020-01-24 [1] standard (@1.1.1)  
##  progress           1.2.2      2019-05-16 [1] CRAN (R 4.1.0)     
##  promises           1.2.0.1    2021-02-11 [1] standard (@1.2.0.1)
##  purrr            * 0.3.4      2020-04-17 [1] standard (@0.3.4)  
##  R6                 2.5.1      2021-08-19 [1] standard (@2.5.1)  
##  RANN               2.6.1      2019-01-08 [1] standard (@2.6.1)  
##  rappdirs           0.3.3      2021-01-31 [1] standard (@0.3.3)  
##  RColorBrewer       1.1-2      2014-12-07 [1] standard (@1.1-2)  
##  Rcpp               1.0.7      2021-07-07 [1] standard (@1.0.7)  
##  RcppAnnoy          0.0.19     2021-07-30 [1] standard (@0.0.19) 
##  RCurl              1.98-1.5   2021-09-17 [1] standard (@1.98-1.)
##  readr            * 2.0.2      2021-09-27 [1] standard (@2.0.2)  
##  readxl             1.3.1      2019-03-13 [1] standard (@1.3.1)  
##  reprex             2.0.1      2021-08-05 [1] standard (@2.0.1)  
##  reshape2           1.4.4      2020-04-09 [1] standard (@1.4.4)  
##  reticulate         1.22       2021-09-17 [1] standard (@1.22)   
##  rlang              0.4.12     2021-10-18 [1] standard (@0.4.12) 
##  rmarkdown          2.11       2021-09-14 [1] standard (@2.11)   
##  ROCR               1.0-11     2020-05-02 [1] standard (@1.0-11) 
##  rpart              4.1-15     2019-04-12 [2] CRAN (R 4.1.0)     
##  rprojroot          2.0.2      2020-11-15 [1] standard (@2.0.2)  
##  RSQLite            2.2.8      2021-08-21 [1] standard (@2.2.8)  
##  rstudioapi         0.13       2020-11-12 [1] standard (@0.13)   
##  Rtsne              0.15       2018-11-10 [1] standard (@0.15)   
##  rvest              1.0.2      2021-10-16 [1] standard (@1.0.2)  
##  S4Vectors          0.32.0     2021-10-26 [1] standard (@0.32.0) 
##  sass               0.4.0      2021-05-12 [1] standard (@0.4.0)  
##  scales             1.1.1      2020-05-11 [1] standard (@1.1.1)  
##  scattermore        0.7        2020-11-24 [1] standard (@0.7)    
##  schelpr          * 0.0.9017   2021-11-30 [1] local              
##  sctransform        0.3.2      2020-12-16 [1] standard (@0.3.2)  
##  sessioninfo        1.1.1      2018-11-05 [2] CRAN (R 4.1.0)     
##  Seurat             4.0.5      2021-10-17 [1] standard (@4.0.5)  
##  SeuratObject       4.0.3      2021-11-10 [1] standard (@4.0.3)  
##  shiny              1.7.1      2021-10-02 [1] standard (@1.7.1)  
##  sir              * 0.0.0.9011 2021-11-26 [1] local              
##  spatstat.core      2.3-1      2021-11-02 [1] standard (@2.3-1)  
##  spatstat.data      2.1-0      2021-03-21 [1] standard (@2.1-0)  
##  spatstat.geom      2.3-0      2021-10-09 [1] standard (@2.3-0)  
##  spatstat.sparse    2.0-0      2021-03-16 [1] standard (@2.0-0)  
##  spatstat.utils     2.2-0      2021-06-14 [1] standard (@2.2-0)  
##  stringi            1.7.5      2021-10-04 [1] standard (@1.7.5)  
##  stringr          * 1.4.0      2019-02-10 [1] standard (@1.4.0)  
##  survival           3.2-13     2021-08-24 [1] standard (@3.2-13) 
##  tensor             1.5        2012-05-05 [1] standard (@1.5)    
##  tibble           * 3.1.6      2021-11-07 [1] standard (@3.1.6)  
##  tidyr            * 1.1.4      2021-09-27 [1] standard (@1.1.4)  
##  tidyselect         1.1.1      2021-04-30 [1] standard (@1.1.1)  
##  tidyverse        * 1.3.1      2021-04-15 [1] standard (@1.3.1)  
##  tweenr             1.0.2      2021-03-23 [1] standard (@1.0.2)  
##  tzdb               0.2.0      2021-10-27 [1] standard (@0.2.0)  
##  utf8               1.2.2      2021-07-24 [1] standard (@1.2.2)  
##  uwot               0.1.10     2020-12-15 [1] standard (@0.1.10) 
##  vctrs              0.3.8      2021-04-29 [1] standard (@0.3.8)  
##  viridisLite        0.4.0      2021-04-13 [1] standard (@0.4.0)  
##  withr              2.4.2      2021-04-18 [1] standard (@2.4.2)  
##  xfun               0.27       2021-10-18 [1] standard (@0.27)   
##  XML                3.99-0.8   2021-09-17 [1] standard (@3.99-0.)
##  xml2               1.3.2      2020-04-23 [1] standard (@1.3.2)  
##  xtable             1.8-4      2019-04-21 [1] standard (@1.8-4)  
##  XVector            0.34.0     2021-10-26 [1] standard (@0.34.0) 
##  yaml               2.2.1      2020-02-01 [1] standard (@2.2.1)  
##  zlibbioc           1.40.0     2021-10-26 [1] standard (@1.40.0) 
##  zoo                1.8-9      2021-03-09 [1] standard (@1.8-9)  
## 
## [1] /home/gleixner/R/x86_64-pc-linux-gnu-library/4.1.0b
## [2] /software/r/4.1.0/lib64/R/library
load helper functions
sgi_blue    = '#5087C8'
sgi_yellow1 = '#F2EE35'
sgi_yellow2 = '#FED98E'
b110_grey   = '#808080'
b110_grey_light   = '#909090'
b110_transparent_black = ggplot2::alpha('#000000', 0.5)
google_red = '#dd4b39'
google_green = '#0F9D58'
google_yellow = '#F4B400'
google_blue = '#4285F4'
b110_colors <- c(sgi_blue, sgi_yellow1, sgi_yellow2, b110_grey, b110_grey_light, b110_transparent_black, google_red, google_green, google_yellow, google_blue)
legend_dot_size <- 1

maja_red <- alpha('#E64B35CC', 0.8)
maja_blue <- alpha('#4DBBD5CC', 0.8)
maja_gray <- alpha('gray', 0.8)

labeller_set <- function(width = Inf, ...) {
  if (width == Inf) {
    labeller(set = labels_of_set, ...) # returned
  } else {
    wrapped_labels_of_set <- stringr::str_wrap(labels_of_set, width = width)
    names(wrapped_labels_of_set)  <- levels_of_set
    labeller(set = wrapped_labels_of_set, ...)  # returned
  }
}

levels_of_cell_type <- c("pseudo bulk", "dead cells", "Stem",       "TA",       "Enterocyte Progenitor", "Enterocyte",  "EEC",     "Paneth",       "Goblet & Paneth",         "Goblet",       "Tuft"      )
labels_of_cell_type <- c("pseudo bulk", "dead cells", "Stem cells", "TA cells", "EC prog.",              "Enterocytes", "EECs",    "Paneth cells", "Goblet &\nPaneth cells",  "Goblet cells", "Tuft cells")
colors_of_cell_type <- c("#303030",     "#303030",    "#66a61e",    "#1b9e77",  "#418DA6",               "#7570b3",     "#a846a0", "#ce5a02",      "#a6761d",                 "#e6ab02",      "#FE5D26"   )  # old tuft: e7298a #E53D00
names(colors_of_cell_type) <- levels_of_cell_type

scale_x_cell_type <- function(
  name = "cell type",
  breaks = levels_of_cell_type,
  labels = labels_of_cell_type,
  limits = function(x) levels_of_cell_type[levels_of_cell_type %in% x],
  guide = guide_axis(
    #n.dodge=2,
    angle = 45,
    title = element_blank()),
  ...
) scale_x_discrete(name = name, breaks = breaks, labels = labels, guide = guide, limits = limits, ...)

scale_y_cell_type <- function(
  name = "cell type",
  breaks = rev(levels_of_cell_type),
  labels = rev(labels_of_cell_type),
  limits = function(x) rev(levels_of_cell_type[levels_of_cell_type %in% x]),
  guide = guide_axis(
    #n.dodge=2,
    #angle = 45,
    title = element_blank()),
  ...
) scale_y_discrete(name = name, breaks = breaks, labels = labels, guide = guide, limits = limits, ...)

scale_color_cell_type <- function(
  name = "cell type",
  breaks = levels_of_cell_type,
  labels = labels_of_cell_type,
  values = colors_of_cell_type,
  limits = identity,
  na.value = "#808080",
  guide = guide_legend(override.aes = list(size = legend_dot_size)),
  ...
) scale_color_manual(
  name = name,
  breaks = breaks,
  limits = limits,
  labels = labels,
  values = values,
  na.value = na.value,
  guide = guide,
  ...
)

scale_fill_cell_type <- function(
  name = "cell type",
  breaks = levels_of_cell_type,
  labels = labels_of_cell_type,
  values = colors_of_cell_type,
  limits = identity,
  na.value = "#808080",
  guide = guide_legend(override.aes = list(size = legend_dot_size)),
  ...
) scale_fill_manual(
  name = name,
  breaks = breaks,
  limits = limits,
  labels = labels,
  values = values,
  na.value = na.value,
  guide = guide,
  ...
)

scale_x_age <- function(..., limits = c("young", "aged")) scale_x_discrete(..., limits = limits)

age_colors <- c(young = "#91D1C2CC", aged = "#3C5488CC")
age_colors_light <- c(young = "#91D1C2CC", aged = "#8491B4CC")
age_colors_dark <- c(young = "#00A087CC", aged = "#3C5488CC")

scale_color_age <- function(..., limits = names(age_colors), values = age_colors, na.value = "grey", guide = guide_legend(override.aes = list(size = legend_dot_size))) ggplot2::scale_color_manual(..., limits = limits, values = values, na.value = na.value, guide = guide)
#scale_color_age_dark <- function(..., values = age_colors_dark) scale_color_age(..., values = values)

scale_fill_age <- function(..., limits = names(age_colors), values = age_colors, na.value = "grey", guide = guide_legend(override.aes = list(size = legend_dot_size))) ggplot2::scale_fill_manual(..., limits = limits, values = values, na.value = na.value, guide = guide)
scale_fill_age_dark <- function(..., values = age_colors_dark) scale_fill_age(..., values = values)

volcano_plot <- function(dt, alpha = 0.1, highlight = rank(padj)<=20, color = "black", box.padding = 0.25, label_size = 4, ...){
  if(missing(alpha)) message(paste0("FDR <= alpha = ", alpha))
  dt <- copy(dt)
  dt[padj > alpha, sir.class := "not significant"]
  dt[padj <= alpha & -log2FoldChange > 0, sir.class := "up regulated"]
  dt[padj <= alpha & -log2FoldChange < 0, sir.class := "down regulated"]
  labelled_subset_dt <- eval(substitute(copy(dt)[!(highlight), gene := ""]))
  
  dt %>% ggplot(aes(x = -log2FoldChange, y = -log10(pvalue))) +
    #geom_vline(xintercept = 0, linetype = "dotted") +
    geom_vline(xintercept = 0, linetype = "dashed") +
    ggrastr::rasterize(geom_point(data = dt[sir.class == "not significant",],aes(color = sir.class), size=0.5)) +
    geom_point(data = dt[sir.class != "not significant",], aes(color = sir.class), size=0.5) +
    scale_color_manual(breaks = c("down regulated", "up regulated", "not significant"), values = c(maja_blue, maja_red, maja_gray)) +
    guides(color = guide_legend(title = element_blank())) +
    labs(x = "log2fc aged/young") +
    # geom_vline(xintercept = c(0.5, -0.5), linetype = "dotted") + geom_hline(yintercept = 2, linetype = "dotted") +
    if(nrow(labelled_subset_dt)>0) ggrepel::geom_label_repel(aes(label = gene), data = labelled_subset_dt, color = color, min.segment.length = 0, size = label_size, box.padding = box.padding, max.overlaps = Inf, ...) else NULL +
    NULL
}

geom_cell <- function(..., size = 0.2, raster.dpi = 600) ggrastr::geom_point_rast(..., size = size, raster.dpi = raster.dpi)
load cell meta data
dt <- fread(fs::path(processed_data_dir, "meta_and_reductions.tsv.gz"), sep= "\t")
tSNE_by_sample_id
dt %>% ggplot(aes(x=tSNE_1, y=tSNE_2, color=sample_id)) +
  geom_cell() +
  # scale_color_sample_id() +
  NULL

tSNE_by_cell_type
dt %>% ggplot(aes(x=tSNE_1, y=tSNE_2, color=cell_type)) +
  geom_cell() +
  scale_color_cell_type() +
  NULL

tSNE_by_age
dt %>% ggplot(aes(x=tSNE_1, y=tSNE_2, color=condition)) +
  geom_cell() +
  scale_color_age() +
  NULL

cell numbers per type
n_cells_dt <- dt[, .(n_cells = .N), by = .(cell_type, sample_id, treatment, condition)][, fraction_of_cells := n_cells/sum(n_cells), by = sample_id]

n_cells_dt[, mean_fraction_of_cells := mean(fraction_of_cells), by= cell_type][, rare := mean_fraction_of_cells < max(mean_fraction_of_cells)/4]

dodge_width <- 0.9
point_jitter_width <- 0.0
this_dodge <- function(width = dodge_width, reverse = TRUE, padding = 0.2) position_dodge2(width = width, reverse = reverse, padding = padding)

p <- n_cells_dt %>% ggplot(aes(x = cell_type, y = fraction_of_cells, color = condition)) +
  stat_summary(aes(fill = condition), fun.data = mean_se, geom = "col", position = this_dodge(), width = dodge_width) +
  stat_summary(aes(color = stage(condition, after_scale = colorspace::darken(color, 0.3))), fun.data = mean_se, geom = "errorbar", position = this_dodge(padding=1/3), width = dodge_width) +
  stat_summary(aes(group = interaction(cell_type, condition), color = stage(condition, after_scale = colorspace::darken(color, 0.3))), fun = mean, geom = "crossbar", position = this_dodge(), width = dodge_width) +
  geom_point(position = position_jitterdodge(point_jitter_width, dodge.width = dodge_width), aes(group = interaction(cell_type, forcats::fct_rev(condition)), fill=condition), color = "black") +
  labs(y = "fraction of cells per sample") +
  guides(color = "none") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_color_age() +
  scale_fill_age() +
  scale_x_cell_type() +
  theme(axis.title.x = element_blank()) +
  NULL
p

cell number proportion change
n_cells_dt[, broom::tidy(t.test(log(fraction_of_cells)~condition, data = .SD)), by = cell_type][, .(relative_change_fraction_ci = pmax(1-exp(conf.low), exp(conf.high)-1)), by = cell_type][]
n_cells_dt[, broom::tidy(t.test(fraction_of_cells~condition, data = .SD)), by = cell_type
           ][, .(mean_fraction_young = estimate2, absolute_change_fraction_ci = pmax(abs(conf.low), abs(conf.high)), conf.low, conf.high, pval = p.value), by = cell_type
             ][, relative_change_fraction_ci := scales::percent_format(1)(absolute_change_fraction_ci/mean_fraction_young)
               ][, adj_pval := p.adjust(pval, method = "BH")][]

DEGs

unnamed-chunk-1
de_single_cell_dt <- fread(fs::path(differential_expression_dir, "DE_samples_as_replicates.tsv.gz"), sep="\t")
de_single_cell_dt <- de_single_cell_dt[contrast == "condition_young_vs_aged"]
de_single_cell_dt[, gene_name := gene]
setnames(de_single_cell_dt, 
         old =  c("estimate", "stderror", "statistic", "p.value", "p.adjusted"),
         new = c("log2FoldChange", "lfcSE", "stat", "pvalue", "padj")
)

cell_types <- de_single_cell_dt[, sort(unique(cell_type))]
cell_types <- c(cell_types[cell_types %in% "pseudo bulk"], cell_types[!(cell_types %in% "pseudo bulk")])

de_single_cell_dt[is.na(padj), padj := 1]
visualize number of degs
n_genes_differentially_expressed_per_celltype <- de_single_cell_dt[cell_type != "pseudo bulk"][padj<0.1][, .N , keyby = cell_type][levels_of_cell_type,][is.na(N), N := 0]
n_genes_differentially_expressed_per_celltype
de_single_cell_dt[cell_type != "pseudo bulk"] %>% ggplot(aes(x = cell_type, y = -log2FoldChange)) +
  ggrastr::rasterize(geom_jitter(color = "gray"), dpi = 600) +
  ggforce::geom_sina(aes(color = cell_type), data = de_single_cell_dt[cell_type != "pseudo bulk"][padj<0.1], scale="width") +
  labs(y = "log2 FC aged / young") +
  guides(color = "none") +
  geom_hline(yintercept = 0, color = "gray50", linetype="dashed") +
  scale_x_cell_type() +
  scale_color_cell_type() +
  theme(axis.title.x = element_blank(), aspect.ratio = 2/3) +
  NULL

# ggsave("results/deg_count.pdf", width = 7.7, height = 3.5, bg = "transparent")
volcano
n_highlight <- 20
volcano_plot(de_single_cell_dt[cell_type=="pseudo bulk",],
  highlight = rank(- (abs(log2FoldChange)-(lfcSE*qnorm(1-0.1*n_highlight/length(lfcSE)))))<=n_highlight, 
  label_size = 2.5, 
  force_pull = 5, 
  point.padding = 0.2,
  force = 70,
  #nudge_y = 10,
  color = "black",     # text color
  bg.color = "white", # shadow color
  bg.r = 0.2          # shadow radius
) +
  facet_wrap("cell_type") +
  NULL
## FDR <= alpha = 0.1
## Warning: Ignoring unknown parameters: bg.colour, bg.r

mhcII pval
mhcII_features_dt <- fread(fs::path(meta_data_dir, "SI_feature_gene_lists.tsv"))[feature %in% "MHC class II"]
mhcII_features_dt <- de_single_cell_dt[mhcII_features_dt, on = .(gene = marker_gene_symbol)]

mhcII_features_dt[cell_type != "pseudo bulk"] %>%
  ggplot(aes(x = gene, y = cell_type, color = -log2FoldChange)) +
  geom_point(aes(size = -log10(padj))) +
  scale_size_area(expression("-log"[10]*"\n p"["adj."]), n.breaks = 3)+
  scale_y_cell_type() +
  scale_color_viridis_c(option = "inferno", direction = -1) +
  guides(color = guide_colorbar(expression("log"[2]*"FC")), x=guide_axis(angle=45)) +
  NULL
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : font metrics unknown for character 0xa

supplementary figure qc filtering
dt %>% ggplot(aes(x=nCount_RNA, y=fraction_mito, color=sample_id)) +
  geom_cell() +
  facet_wrap("cell_type", nrow=2) +
  labs(x = "counts", y = "fraction of counts from mitochondrial genes") +
  theme(
    panel.grid.major.x = element_line(color="gray95", size = 0.3),
    panel.grid.minor.x = element_line(color="gray95", size = 0.1),
    panel.grid.major.y = element_line(color="gray95", size = 0.3),
    panel.grid.minor.y = element_line(color="gray95", size = 0.1)
  ) +
  schelpr::scale_x_log()+
  guides(color = guide_legend(override.aes = list(size = 1))) +
  NULL

inflammaging score tsne
dt[, age := factor(condition, c("young", "aged"))]

p1 <- dt %>% ggplot(aes(x=tSNE_1, y=tSNE_2, color=inflammaging_score_scaled)) + 
  geom_cell() + 
  labs(color="inflammaging score") +
  NULL
p1 + facet_wrap("age") + scale_color_gradientn(colours = heatmaply::cool_warm(500), limits = c(0, 1))

dt %>% ggplot(aes(x=tSNE_1, y=tSNE_2, color=rank(inflammaging_score_scaled)/length(inflammaging_score_scaled))) + 
  geom_cell() + 
  labs(color="inflammaging score\n(quantile normalized)") +
  NULL + facet_wrap("age") + scale_color_gradientn(colors = heatmaply::cool_warm(500), limits= c(0,1), labels=scales::percent_format())

inflammaging score
dts <- dt[, .(inflammaging_score = mean(inflammaging_score), inflammaging_score_scaled = mean(inflammaging_score_scaled)), by = .(cell_type, condition, Tx_run_number, sample_id)]
dts_wide <- dcast(dts, cell_type + Tx_run_number ~ condition, value.var = c("inflammaging_score", "inflammaging_score_scaled"))
dts_wides <- dts_wide[, do.call(cbind, lapply(.SD, mean_cl_normal)), .SDcols = c("inflammaging_score_scaled_young", "inflammaging_score_scaled_aged"),  by= cell_type]

dts %>%
  ggplot(aes(x = inflammaging_score_scaled, y = cell_type, color = condition)) +
  xlab("inflammaging score") +
  #geom_line(aes(group = cell_type), color = "gray", size = 1) +
  #stat_summary(aes(group = stage(interaction(cell_type, condition), after_stat = cell_type)), geom = "line") +
  geom_segment(data = dts_wides, aes(x = inflammaging_score_scaled_young.y, xend = inflammaging_score_scaled_aged.y, yend = cell_type), color = "gray",  size = 1) +
  stat_summary(geom = "point", size = 5) +
  geom_point() +
  #geom_point(aes(x = cell_type_inflammaging_score), size = 2, alpha = 0.5) +
  #cowplot::theme_cowplot() +
  scale_color_age() +
  scale_y_cell_type() +
  NULL
## No summary function supplied, defaulting to `mean_se()`